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ABSTRACT 

A global nonlinear parametric modeling technique is 
described and demonstrated. The technique uses multivariate 
orthogonal modeling functions generated from the data to 
determine nonlinear model structure, then expands each 
retained modeling function into an ordinary multivariate 
polynomial. The final model form is a finite multivariate 
power series expansion for the dependent variable in terms 
of the independent variables. Partial derivatives of the 
identified models can be used to assemble globally valid 
linear parameter varying models. The technique is 
demonstrated by identifying global nonlinear parametric 
models for nondimensional aerodynamic force and moment 
coefficients from a subsonic wind tunnel database for the 
F-16 fighter aircraft. Results show less than 10% difference 
between wind tunnel aerodynamic data and the nonlinear 
parameterized model for a simulated doublet maneuver at 
moderate angle of attack. Analysis indicated that the global 
nonlinear parametric models adequately captured the 
multivariate nonlinear aerodynamic functional dependence. 

NOMENCLATURE 

wing reference chord, ft 
wing span, ft 

aerodynamic force coefficients 
aerodynamic moment coefficients 
body axis angular rates, rad/sec 
predicted squared error 
airspeed, ft/sec 
angle of attack, rad 
sideslip angle, rad 

elevator, aileron, rudder deflections, rad 
longitudinal c.g. position 
reference longitudinal c.g. position = 0.35 

1. INTRODUCTION 

An important aspect of accurately modeling nonlinear 
functional dependence is determining the mathematical form 
relating independent variables to a dependent variable. In 
general, the goal is to find a compact model structure which 
still has adequate complexity to capture the nonlinearities. 
Keeping the number of terms in the model low improves 
model parameter identifiability, resulting in a more accurate 
model with good prediction capability. 

Models can be loosely classified as local or global. 
Local models are identified using data from a relatively 


small region of the independent variable space. It follows 
that local models are valid for small ranges of the 
independent variables. A global model results when the 
range of validity for the identified model covers a large 
portion of the independent variable space. Given this latter 
type of (generally nonlinear) model, operations can often be 
streamlined by replacing many local models with a single 
global model. 

When the global nonlinear model can be identified in a 
parametric form using simple analytic terms, it is possible 
to formulate Linear Parameter Varying (LPV) models that 
are globally valid. Partial derivatives of the global 
analytical models with respect to the independent variables 
provide global linear model parameter variations for LPV 
models. Such models are useful in robust and nonlinear 
control design. This type of model also provides insight 
not available from a family of linear constant coefficient 
models obtained from local linearization of the nonlinear 
functional dependence at various operating points. 
Calculation of finite difference linear constant coefficient 
models requires that perturbation size be chosen carefully to 
ensure proper linear characterization of a nonlinear 
functional dependency and to avoid inaccuracies due to local 
measurement noise. These considerations are avoided with a 
global nonlinear model. Globally valid analytical models 
and their associated smooth gradients are also useful for 
optimization and global nonlinear stability and control 
analysis. 

Recently, a method has been developed which identifies 
an adequate multivariate polynomial model structure with 
accurate parameter estimates based on orthogonal modeling 
functions generated from the data ' . Selected orthogonal 
modeling functions are included in the model based on 
minimizing the predicted squared error, and then are 
decomposed into ordinary multivariate polynomials. The 
final identified model consists of selected terms from a 
multivariable power series expansion for the dependent 
variable in terms of the independent variables. 

The purpose of this work was to investigate the 
suitability of the orthogonal function modeling technique 
for global nonlinear modeling on a realistic problem, 
namely a nonlinear aerodynamic database- obtained from 
wind tunnel tests of a modern fighter^. The investigation 
focuses on the degree to which a realistic nonlinear 
functional dependence embodied in large tables of values can 
be modeled using compact multivariate polynomial 
expressions while retaining good predictive capability. The 
next section outlines the theoretical development. 

Following this, the nonlinear modeling technique is applied 
to a wind tunnel aerodynamic database for the F-16 aircraft. 
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2. THEORETICAL DEVELOPMENT 

Assume an /V-dimensional vector of dependent variable 
T 

values, j = [yj, J'iv] > modeled in terms of a linear 
combination of n modeling functions pj, j = 1 , 2 ,...,;;. 
Each p j is an ;V-dimensional vector which in general 
depends on the independent variables. Then, 

y = c 1 p l +c 2 p 2 +- + c„p n +z (1) 

The c j , j = 1,2,...,//, are constant model parameters and e 

denotes the modeling error vector. We put aside for the 
moment the question of how to determine the modeling 
functions p j , as well as how to select which functions 

should be included in the model of Eq. (1) (which implicitly 
determines n). Now define an Nxn matrix P, 

P = [P\’Pl’-’Pn\ ( 2 ) 

T 

and let c = Cj,---, c^y] . Then using Eq. (1), 

y = P c + e ( 3 ) 

The goal is to determine c which minimizes the least 
squares cost function 

J = (y-Pc) T {y-Pc) (4) 

The least squares estimate of c is computed from^ 

c = [p 7 'p]” 1 P T y (5) 

and the estimated parameter covariance matrix is^ 

var(c A )=£{(c-c)(c-c) r } = -^^[p r p] 1 (6) 

where E is the expectation operator, and J is the cost 
calculated from Eq. (4) with c = c . 

Reference [1] describes a procedure for using the 
independent variable data to generate orthogonal modeling 
functions, which have the following important property: 

pJPj = 0 > i*j > i,j = l%...,n (7) 

Using Eqs. (2) and (7) in Eq. (5), the jth element of the 
estimated parameter vector c is given by 



Combining Eqs. (2), (4), (7), and (8), 

J=y T y-t(pjyf/(pJpj) ( 9 ) 

7=1 

Eq. (9) shows that when the modeling functions are 
orthogonal, the reduction in the estimated cost resulting 


from including the term c ■ p • in the model depends only 

the dependent variable datay and the added orthogonal 
modeling function pj . This decouples the least squares 

estimation, and makes it possible to evaluate each 
orthogonal modeling function in terms of its ability to 
reduce the least squares model fit to the data, regardless of 
which other orthogonal modeling functions are present in 
the model. The orthogonal modeling functions are chosen 
to minimize predicted squared error PSE , defined by-* 

PSE = — + a 2 — (10) 

N ° N 

where 

1 N r. 1 N 

a l = ^i{yi-y) ’ y = ^i >< (11 > 

v /= 1 v i = 1 

In Eq. (10), the PSE depends on the mean square fit 
error j/N, and a term proportional to the number of terms 
in the model, n . The latter term prevents overfitting with 
too many model terms, which is detrimental to model 
prediction accuracy^. Note that while the mean square fit 
error J/N must decrease with the addition of each 
orthogonal modeling function by Eq. (9), the oveifit penalty 
term erf n/N increases with each added model term (n 

increases), so that PSE always has a single global 
minimum value. Ref. [5] contains details on the statistical 
properties of the PSE metric, including justification for its 
use in modeling problems. 

The orthogonal functions are generated in a manner that 
allows them to be decomposed without ambiguity into an 
expansion of ordinary multivariate polynomials^. The 
process can be repeated to generate orthogonal functions of 
arbitrary order in the independent variables, subject only to 
limitations related to the information contained in the data. 

Using orthogonal functions to model the dependent 
variable made it possible to evaluate the merit of including 
each modeling function individually as part of the model, 
using the predicted squared error, PSE . This approach made 
model structure determination a well-defined and 
straightforward process. After the orthogonal modeling 
functions that minimized PSE were selected, each retained 
orthogonal function was expanded into an ordinary 
polynomial expression, and common terms in the ordinary 
polynomials were combined using double precision 
arithmetic to arrive finally at a multivariate model using 
only ordinary polynomials in the independent variables. 
Ordinary polynomial coefficients with absolute value less 
than 10"^ were dropped from the final model. 

Orthogonal modeling functions are useful in 
determining the model structure for the dependent variable 
using the PSE metric, by virtue of the benefits of 
orthogonal functions and the resultant decoupling of the 
associated least squares problem. The subsequent 
decomposition of the retained orthogonal functions is done 
to express the results in physically meaningful terms and to 
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allow analytic differentiation for partial derivatives of the 
dependent variable with respect to the independent variables. 

3. RESULTS 

Wind tunnel aerodynamic data for a 16% scale model of 
the F-16 aircraft flying at relatively low Mach numbers 
(< 0.6), out of ground effect, with landing gear retracted and 
no external stores, is given in Reference [3]. The data used 
in this work was a slightly simplified version- of the 
original wind tunnel database. Nondimensional coefficients 
that vary nonlinearly with flow angles ( a, /j ) , aircraft 


the resulting identified model structure, with estimated 
parameter values from Table 3. 

Model structures for all the component functions in 
Eqs. (12)-(17) appear in Table 2. Estimated parameter 
values are given in Table 3. Global LPV models can be 
assembled via analytic partial differentiation of these 
polynomial aerodynamic models in combination with the 
equations of motion. 

Figure 1 shows results from a simulated 
lateral/directional doublet maneuver at 10° angle of attack, 

Mach 0.26 at sea level, with x =0.25. The control 

C 6 


angular velocities ( p , q, r), and control surface deflections 
(S e , 8 a , S f . ) characterize the aerodynamic forces and 
moments acting on the aircraft. Dependence of the 
nondimensional coefficients on a is included in the q 
dependencies, due to the manner in which the data is 
collected in the wind tunnel. Each nondimensional 
aerodynamic force and moment coefficient was built up 
from a set of component functions, where each component 
function was determined by a table look-up in the wind 
tunnel database. The expressions for the nondimensional 
aerodynamic force and moment coefficients were 

C x = C x {a,S e ) + C (a)q ( 12 ) 

C 1 

C y = C y (p,8 a ,8 r ) +C yp (a)p+C yr (a)F ( 13 ) 

C z = C z (a,p,8 e ) + C z (a)q ( 14 ) 

Q 

C t = C t (a,p) + C t {a)p + C tr (a)r 

+ C f§ (a,p)8 +C fs (a,p)8 ' L '’ 

u a r 

C m = C m{ a ’ S e) +C m q ^ 8 + C z[ x cg ref ~ x cg] (16 ) 
C n= C n( a ’P) + C n p ( a )P + C n } . ( ®) ? 

+ C (a, ft) 8 + C n g (a, ft) 8 ( 17 ) 

U a u r 


surface inputs shown in the first two plots of Figure 1 were 
taken from actual flight test data for the F-16 aircraft. The 
lower plots in Figure 1 show time histories for sideslip 
angle, roll rate, and nondimensional yawing moment 
coefficient, obtained by applying the control inputs shown 
to a full nonlinear F-16 simulation- using the wind tunnel 
aerodynamic database (solid lines) and the identified 
polynomial aerodynamic models (dashed lines). Fess than 
10% difference was seen in these time histories, indicating 
that the polynomial modeling was successful in capturing 
the nonlinear aerodynamic functional dependence. Time 
histories for the other aircraft responses and nondimensional 
coefficients exhibited a similar level of agreement. 



C y( Xc Sref Xc 8){b) 

where 

p = pb/2V q = qc/2V r = rb/2V (18) 

Each function in Eqs. (12)-( 1 7) was modeled 
individually using the orthogonal function modeling 
technique described above for the independent variable ranges 
shown in Table 1. These independent variable ranges 
represent all the data available in the wind tunnel database. 

For example, the C x (a, 8 g ) function was modeled 



using tabulated values of that function as the dependent 
variable and corresponding tabulated values of a and 8 g as 
the independent variables. The first row of Table 2 specifies 


0 2 4 6 8 10 

Time (sec) 

Figure 1 Lateral Directional Doublet Maneuver 
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4. CONCLUDING REMARKS 

Multivariate orthogonal functions generated from the 
data were used to construct global analytical models for 
nondimensional aerodynamic force and moment coefficients 
of the F-16 aircraft based on a subsonic wind tunnel 
database. Each model was a single ordinary multivariate 
polynomial in the independent variables, valid for the entire 
flight envelope encompassed by the wind tunnel data. For a 
realistic simulated lateral/directional doublet maneuver, 
aircraft response variables and aerodynamic coefficients 
computed using the identified polynomial models matched 
those obtained using the wind tunnel database within 10 %. 
Global nonlinear aerodynamic models like those identified 
here are useful in many applications, including flight 
simulation, control system design, and dynamic analysis. 

The modeling technique described and demonstrated here 
is general and can be applied to data from other physical 
systems. The final result is a compact, global analytical 
model of the nonlinear functional dependence embodied in 
the data, with good predictive capabilities. Smooth global 
analytic derivatives of any order can be easily calculated. 
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Table 1 - Independent Variable Ranges for a Compact 
Global Nonlinear Aerodynamic Model of the F-16 


Lower Bound 

Variable 

Lower Bound 

-0.1745 rad 

a 

0.7854 rad 

(-10 deg) 

(45 deg) 

-0.5236 rad 

P 

0.5236 rad 

(-30 deg) 

(30 deg) 

-0.4363 rad 

8 

0.4363 rad 

(-25 deg) 

v e 

(25 deg) 

-0.3752 rad 

8 

0.3752 rad 

(-21.5 deg) 

^ a 

(21.5 deg) 

-0.5236 rad 

8 

0.5236 rad 

(-30 deg) 

r 

(30 deg) 


Table 2 - Model Structure for a Compact Global 


Nonlinear Aerodynamic Model of the F-16 
Function Model Structure 


C x {a,S e ) 

Gq + ci^OC + a 3 ^e 

+a^oc 8 e + a^a 2 + a^a 2 

C x (a) 
q 

+ b^cc + + b^oc 2 + b^oc^ 

Cy{P,S a ,S r ) 

c oP + c l <>a + c 2 *7 

c y p (a) 

<Jq + <7j a + d 2 a 2 + d-y a 3 

Cy r (a) 

+ 6 ^ OC + £3 

C z (a,p.S e ) 

( /o . +/ 5“7 2 4 2 ] ('-PYh*' 

l+/3« +/ 4 « J 

C z («) 

q 

g 0 +Si« + «2 a 2 + .?3“ 3 + «4 “ 4 

Ci(a.P ) 

/zq P + oc(5 + /?2 oc 2 P + / 13 /J 2 
+h 4 ap 2 + h 5 a 3 p + h 6 oc^p + hj a 2 p 2 

C, (a) 
P 

/q + /j OC + z*2 (X 2 “1” *3 (X 2 

C, (a) 

l r 

j() + j\ (X + J 2 «“ + J 3 « 3 + 74 « 4 

Cu («, P) 

°a 

&0 + cc + k <2 P k ^ cc 2 
+k 4 ap + k 5 a 2 P + k 6 a 2 

Ci (■ a,p ) 

°r 

l^ + l^ot + I 2 P + I 3 O 1 P 
+ 0 c 2 p + cc 2 P + IfoP 2 

C m {a, 8 e ) 

niQ + mj a + m 2 8 e + m 3 a 8 e + m 4 <5 2 
+mj a 2 S e + rn.fi <5 3 + m-y a S 2 

C mq (a) 

Uq + ?q a + «2 a 2 + «3 a 2 
+n 4 a 4 + tifi a 2 

C n (a,P) 

oqP + o^ a/3 + o 2 p 2 + 03 ap 2 
+ o 4 a 2 p + o^a 2 p 2 + o 6 a 3 p 

C n («) 

rip 

Pq+ p\ a + pj a 2 + p-y « 3 + P 4 « 4 


< 7 q + a + q~, a~ 

{a,P) 

S a 

r^+ r^a + r~, p + r^ap 
+ r 4 a 2 p + r 5 a 2 p + r 6 a 2 
+ r 7 a 3 + rg/3 3 + r 9 ap 3 

c n „ (a.p) 

S r 

5q+ oc + ^2 P s 2 ocP 

+ s ^ oc 2 P + cc 2 
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Table 3 - Parameter Values for Global Nonlinear F-16 Aerodynamic Model 


a 0 

-1.943367e-02 

a \ 

2.136104e-01 

a 2 

-2.903457e-01 

a 3 

-3.348641e-03 

a 4 

-2.060504e-01 

a 5 

6.988016e-01 

a 6 

-9.035381e-01 

K 

4.833383e-01 

A 

8.644627e+00 

b 2 

1.131098e+01 

A 

-7.422961e+01 

K 

6.075776e+01 

c 0 

-1.145916e+00 

Cl 

6.016057e-02 

C 2 

1.642479e-01 

d() 

-1.006733e-01 

d x 

8.679799e-01 

d 2 

4.260586e+00 

d 3 

-6.923267e+00 

e o 

8.071648e-01 

Cl 

1.189633e-01 

c 2 

4.177702e+00 

c 3 

-9.162236e+00 


-1.378278e-01 

/, 

^1.211369e+00 

A 

4.775 187e+00 

A 

-1.026225e+01 

A 

8.399763e+00 

A 

^f.354000e-01 

So 

-3.054956e+01 

81 


82 

3.292788e+02 

83 

-6.848038e+02 

84 

4.080244e+02 


K 

-1.058583e-01 

K 

-5.776677e-01 

h 2 

-1.672435e-02 

h 3 

1.357256e-01 


2.172952e-01 

mm 

3.464156e+00 


-2.83545 le+00 

A 

-1.098104e+00 

h 

^1.126806e-01 

h 

-1.189974e-01 

h 

1.24772 le+00 

h 

-7.391 132e-01 

Jo 

6.250437e-02 

h 

6.067723e-01 

Ji 

-1.101964e+00 

h 

9.100087e+00 


-1.192672e+01 


-1.463144e-01 

K 

^1.073901e-02 

k 2 

3.253 159e-02 

k 3 

4.851209e-01 

k 4 

2.978850e-01 


-3.746393e-01 

K 

-3.213068e-01 

A 

2.635729e-02 

k 

-2.192910e-02 

i 2 

-3.152901e-03 

h 

-5.817803e-02 

h 

4.516159e-01 

h 

^1.928702e-01 

h 

-1.579864e-02 

m 0 

-2.029370e-02 

m 1 

4.660702e-02 

m 2 

-6.012308e-01 

m 3 

-8.062977e-02 

m 4 

8.320429e-02 

m 5 

5.018538e-01 

m 6 

6.378864e-01 

m 7 

4.226356e-01 


n o 

-5.159153e+00 

«i 

-3.554716e+00 

n 2 

-3.598636e+01 

n 3 

2.247355e+02 

n 4 

^1.1 2099 le+02 

n 5 

2.41 1750e+02 

°0 

2.993363e-01 


6.594004e-02 

o 2 

-2.003 125e-01 

°3 

-6.233977e-02 

O 4 

-2.107885e+00 

°5 

2.141420e+00 

°6 

8.476901e-01 

P() 

2.677652e-02 

Pi 

-3.298246e-01 


1.926178e-01 


4.013325e+00 


^1.404302e+00 


-3.698756e-01 


-1.167551e-01 


-7.641297e-01 


-3.348717e-02 

r i 

4.276655e-02 

— i 

6.573646e-03 

mm 

3.535831e-01 

r 4 

-1.373308e+00 


1.237582e+00 


2.302543e-01 

r i 

-2.512876e-01 

mm 

1.588105-01 

mm 

-5.199526e-01 

s 0 

-8.115894e-02 

s i 

-1.156580e-02 

* 2 

2.514167e-02 

5 3 

2.038748e-01 

EMU 

-3.337476e-01 

5 5 

1.004297e-01 
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